/*
 * tJU-SGA.cpp
 *
 *  Created on: 16-07-2012
 *      Author: mith
 */

#include<cstdlib>
#include<cstdio>
#include<iostream>
#include<math.h>
#include<cstring>
#include<time.h>
#include<gsl/gsl_vector.h> // GSL LIBRARY - remember use flag -lgsl -lgslcblas
#include<gsl/gsl_multiroots.h>
#include<gsl/gsl_errno.h> // to interpolation
#include<gsl/gsl_spline.h>
#include<omp.h> // flaga -fopenmp
#include"tJU-SGA.h" // function definitions here!
using namespace std;

// DEFINITIONS
#define MAX_ITERATION 200
#define PRECISION 1e-5
#define NO_OF_EQUATIONS 8
#define DEBUG 1
#define DEEP_DEBUG 0

/*
 *********************************************************
 ***************** CLASS Functions ***********************
 *********************************************************
*/

class Func {
private:
	parameters p; // parameters  p.n  p.J  p.t  p.tPrim  p.U  p.beta  p.L
	variables v; // variables  v.chiAB  v.chiS  v.chiT  v.DeltaS  v.DeltaT  v.mu  v.mAF  v.d
	double n, r, w, d, m, mu; // contractions

	// Abbreviations - dependend on kx, ky
	double ek;
	double ekPrim;
	double etak;

	double K1;
	double dK1dlchiAB;
	double dK1dlchiS;
	double dK1dlchiT;
	double dK1dlDeltaS;
	double dK1dlDeltaT;
	double dK1dln;
	double dK1dlmAF;

	double K2;
	double dK2dlchiAB;
	double dK2dlchiS;
	double dK2dlchiT;
	double dK2dlDeltaS;
	double dK2dlDeltaT;
	double dK2dln;
	double dK2dlmAF;

public:
	double E[4];
	double dEdlchiAB[4];
	double dEdlchiS[4];
	double dEdlchiT[4];
	double dEdlDeltaS[4];
	double dEdlDeltaT[4];
	double dEdln[4];
	double dEdlmAF[4];

	double gt, gs, dgtdm, dgtdd, dgtdn, dgsdm, dgsdd, dgsdn;
	double lchiAB, lchiS, lchiT, lDeltaS, lDeltaT, ln, lmAF;

	Func(const parameters &p, const variables &v) {
		set(p, v);
		refresh();
	}

	void set(const parameters &p, const variables &v) {
		copyTo(this->p, p);
		copyTo(this->v, v);
		refresh();
	}

	void set(const parameters &p) {
		copyTo(this->p, p);
		refresh();
	}

	// Gutzwiller
	double fgt();
	double fgs();
	double fdgtdm();
	double fdgtdd();
	double fdgtdn();
	double fdgsdm();
	double fdgsdd();
	double fdgsdn();

	// DEBUG
	void printLambda() {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e (l)\n", lchiAB, lchiS, lchiT, lDeltaS, lDeltaT, ln, lmAF);
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e (g)\n", gt, gs, dgtdm, dgtdd, dgtdn, dgsdm, dgsdd, dgsdn);
	}

	// REFRESH
	void refresh() {
		n = p.n;
		r = p.n/2. + v.mAF;
		w = p.n/2. - v.mAF;
		d = v.d;
		mu = v.mu;
		m = v.mAF;

		// obliczam funkcje Gutzwillara i poszczegolne pochodne
		gt = fgt();
		gs = fgs();
		dgtdn = fdgtdn();
		dgtdd = fdgtdd();
		dgtdm = fdgtdm();
		dgsdn = fdgsdn();
		dgsdm = fdgsdm();
		dgsdd = fdgsdd();

		// obliczanie poszczegolnych lambda:
		lchiAB = -gt*p.t + (3./4.)*gs*p.J*v.chiAB;
		lchiS = -gt*p.tPrim;
		lchiT = 0.;
		lDeltaS = -(1./4.)*gs*p.J*v.DeltaS;
		lDeltaT = -(3./4.)*gs*p.J*v.DeltaT;
		ln = -16.*(p.t*v.chiAB + p.tPrim*v.chiS)*dgtdn - \
			2.*p.J*(-3.*v.chiAB*v.chiAB + v.DeltaS*v.DeltaS + 3.*v.DeltaT*v.DeltaT + v.mAF*v.mAF/2.)*dgsdn;
		lmAF = -16.*(p.t*v.chiAB + p.tPrim*v.chiS)*dgtdm - \
			2.*p.J*(-3.*v.chiAB*v.chiAB + v.DeltaS*v.DeltaS + 3.*v.DeltaT*v.DeltaT + v.mAF*v.mAF/2.)*dgsdm - \
			2.*v.mAF*gs*p.J;

	}

	//REFRESH XY
	void refresh(double kx, double ky) {
		ek = 2.*(cos(kx) + cos(ky));
		ekPrim = 4.*cos(kx)*cos(ky);
		etak = 2.*(cos(kx) - cos(ky));

		refresh();
		K1 = 4.*Power(ek,2.)*Power(lchiAB,2.) + 4.*Power(ekPrim,2.)*(Power(lchiS,2.) + Power(lchiT,2.)) + 4.*Power(etak,2.)*Power(lDeltaS,2.) + 4.*Power(etak,2.)*Power(lDeltaT,2.) + Power(lmAF,2.) + Power(ln,2.) + 4.*ln*mu + 4.*Power(mu,2.) + 4.*ekPrim*(lchiT*lmAF + lchiS*ln + 2.*lchiS*mu);
		dK1dlchiAB = 8.*Power(ek,2.)*lchiAB;
		dK1dlchiS = 4.*ekPrim*(2.*ekPrim*lchiS + ln + 2.*mu);
		dK1dlchiT = 4.*ekPrim*(2.*ekPrim*lchiT + lmAF);
		dK1dlDeltaS = 8.*Power(etak,2.)*lDeltaS;
		dK1dlDeltaT = 8.*Power(etak,2.)*lDeltaT;
		dK1dln = 2.*(2.*ekPrim*lchiS + ln + 2.*mu);
		dK1dlmAF = 2.*(2.*ekPrim*lchiT + lmAF);

		K2 = 4.*Power(ek,2.)*Power(lchiAB,2.)*(4.*Power(ekPrim,2.)*Power(lchiS,2.) + 4.*Power(etak,2.)*Power(lDeltaT,2.) + 4.*ekPrim*lchiS*(ln + 2.*mu) + Power(ln + 2.*mu,2.)) + Power(4.*Power(ekPrim,2.)*lchiS*lchiT + 4.*Power(etak,2.)*lDeltaS*lDeltaT + lmAF*(ln + 2.*mu) + 2.*ekPrim*(lchiS*lmAF + lchiT*(ln + 2.*mu)),2.);
		dK2dlchiAB = 8.*Power(ek,2.)*lchiAB*(4.*Power(ekPrim,2.)*Power(lchiS,2.) + 4.*Power(etak,2.)*Power(lDeltaT,2.) + 4.*ekPrim*lchiS*(ln + 2.*mu) + Power(ln + 2.*mu,2.));
		dK2dlchiS = 4.*ekPrim*(4.*Power(ek,2.)*Power(lchiAB,2.)*(2.*ekPrim*lchiS + ln + 2.*mu) + (2.*ekPrim*lchiT + lmAF)*(4.*Power(ekPrim,2.)*lchiS*lchiT + 4.*Power(etak,2.)*lDeltaS*lDeltaT + lmAF*ln + 2.*lmAF*mu + 2.*ekPrim*(lchiS*lmAF + lchiT*ln + 2.*lchiT*mu)));
		dK2dlchiT = 4.*ekPrim*(2.*ekPrim*lchiS + ln + 2.*mu)*(4.*Power(ekPrim,2.)*lchiS*lchiT + 4.*Power(etak,2.)*lDeltaS*lDeltaT + lmAF*(ln + 2.*mu) + 2.*ekPrim*(lchiS*lmAF + lchiT*(ln + 2.*mu)));
		dK2dlDeltaS = 8.*Power(etak,2.)*lDeltaT*(4.*Power(ekPrim,2.)*lchiS*lchiT + 4.*Power(etak,2.)*lDeltaS*lDeltaT + lmAF*(ln + 2.*mu) + 2.*ekPrim*(lchiS*lmAF + lchiT*(ln + 2.*mu)));
		dK2dlDeltaT = 8.*Power(etak,2.)*(4.*Power(ek,2.)*Power(lchiAB,2.)*lDeltaT + lDeltaS*(4.*Power(ekPrim,2.)*lchiS*lchiT + 4.*Power(etak,2.)*lDeltaS*lDeltaT + lmAF*(ln + 2.*mu) + 2.*ekPrim*(lchiS*lmAF + lchiT*(ln + 2.*mu))));
		dK2dln = 8.*Power(ek,2.)*Power(lchiAB,2.)*(2.*ekPrim*lchiS + ln + 2.*mu) + 2.*(2.*ekPrim*lchiT + lmAF)*(4.*Power(ekPrim,2.)*lchiS*lchiT + 4.*Power(etak,2.)*lDeltaS*lDeltaT + lmAF*ln + 2.*lmAF*mu + 2.*ekPrim*(lchiS*lmAF + lchiT*ln + 2.*lchiT*mu));
		dK2dlmAF = 2.*(2.*ekPrim*lchiS + ln + 2.*mu)*(4.*Power(ekPrim,2.)*lchiS*lchiT + 4.*Power(etak,2.)*lDeltaS*lDeltaT + lmAF*(ln + 2.*mu) + 2.*ekPrim*(lchiS*lmAF + lchiT*(ln + 2.*mu)));

		E[0]          = -0.5*sqrt(K1 - 2.*sqrt(K2));
		dEdlchiAB[0]  = (dK1dlchiAB - dK2dlchiAB/sqrt(K2))/(8.*E[0]);
		dEdlchiS[0]   = (dK1dlchiS - dK2dlchiS/sqrt(K2))/(8.*E[0]);
		dEdlchiT[0]   = (dK1dlchiT - dK2dlchiT/sqrt(K2))/(8.*E[0]);
		dEdlDeltaS[0] = (dK1dlDeltaS - dK2dlDeltaS/sqrt(K2))/(8.*E[0]);
		dEdlDeltaT[0] = (dK1dlDeltaT - dK2dlDeltaT/sqrt(K2))/(8.*E[0]);
		dEdln[0]      = (dK1dln - dK2dln/sqrt(K2))/(8.*E[0]);
		dEdlmAF[0]    = (dK1dlmAF - dK2dlmAF/sqrt(K2))/(8.*E[0]);

		E[1]          = -E[0];
		dEdlchiAB[1]  = -dEdlchiAB[0];
		dEdlchiS[1]   = -dEdlchiS[0];
		dEdlchiT[1]   = -dEdlchiT[0];
		dEdlDeltaS[1] = -dEdlDeltaS[0];
		dEdlDeltaT[1] = -dEdlDeltaT[0];
		dEdln[1]      = -dEdln[0];
		dEdlmAF[1]    = -dEdlmAF[0];

		E[2]          = -0.5*sqrt(K1 + 2.*sqrt(K2));
		dEdlchiAB[2]  = (dK1dlchiAB + dK2dlchiAB/sqrt(K2))/(8.*E[2]);
		dEdlchiS[2]   = (dK1dlchiS + dK2dlchiS/sqrt(K2))/(8.*E[2]);
		dEdlchiT[2]   = (dK1dlchiT + dK2dlchiT/sqrt(K2))/(8.*E[2]);
		dEdlDeltaS[2] = (dK1dlDeltaS + dK2dlDeltaS/sqrt(K2))/(8.*E[2]);
		dEdlDeltaT[2] = (dK1dlDeltaT + dK2dlDeltaT/sqrt(K2))/(8.*E[2]);
		dEdln[2]      = (dK1dln + dK2dln/sqrt(K2))/(8.*E[2]);
		dEdlmAF[2]    = (dK1dlmAF + dK2dlmAF/sqrt(K2))/(8.*E[2]);

		E[3]          = -E[2];
		dEdlchiAB[3]  = -dEdlchiAB[2];
		dEdlchiS[3]   = -dEdlchiS[2];
		dEdlchiT[3]   = -dEdlchiT[2];
		dEdlDeltaS[3] = -dEdlDeltaS[2];
		dEdlDeltaT[3] = -dEdlDeltaT[2];
		dEdln[3]      = -dEdln[2];
		dEdlmAF[3]    = -dEdlmAF[2];
	}

   void onePointE(double w, int i, int j) {
      double x = (double)i*2.*M_PI/p.L;
      double y = (double)j*2.*M_PI/p.L;
      refresh(x, y);
	      for(int k=0; k<4; k++) {
         double ebE = exp(-p.beta*E[i]);
         if(ebE > 100000) { w += -p.beta*E[i]; }
         else if(ebE < 0.001) { w += ebE; }
         else { w += log(1. + exp(-p.beta*E[i])); }
      }
   }

   double Energy() {
      //refresh();
      double en = 0;
      double temp = 0;
      unsigned int max = p.L/2; // attention: L should be even!

      // SERIAL CASE
      // SMART WAY OF SUMMATION - I use symmetry k_x,y <-> -k_x,y
      // * edges (2 times)
      for(unsigned int i=1; i<max; i++) {
         onePointE(temp, i, 0);
         onePointE(temp, 0, i);
         onePointE(temp, i, max-i);
      }
      en += temp*2.; temp=0;
      // * interior area (4 times) - SERIAL CASE
      for(unsigned int i=1; i<max; i++) {
         for(unsigned int j=1; j<i; j++) {
            onePointE(temp, max-i, j);
         }
      }
      en += temp*4.; temp=0;
      // * missing corners
      onePointE(en, max, 0);
      onePointE(en, 0, max);
      onePointE(en, 0, 0);

      refresh();
      en /= -p.beta*p.L*p.L;
      en += 8.*gt*(p.t*v.chiAB + p.tPrim*v.chiS) + p.U*v.d*v.d \
         + gs*p.J*(-3.*v.chiAB*v.chiAB + v.DeltaS*v.DeltaS + 3.*v.DeltaT*v.DeltaT + v.mAF*v.mAF/2.) \
         + 0.5*(ln*p.n + lmAF*v.mAF) \
         + 8.*(lchiAB*v.chiAB + lchiS*v.chiS + lchiT*v.chiT + lDeltaS*v.DeltaS + lDeltaT*v.DeltaT) \
         + v.mu*p.n;
      return en;
   }


};

double Func::fgt() {
	double part1 = (n-2.*v.d)/(n-2.*r*w);
	double part2 = sqrt((1.-w)*(1.-n+v.d)/(1.-r)) + sqrt(w*v.d/r);
	double part3 = sqrt((1.-r)*(1.-n+v.d)/(1.-w)) + sqrt(r*v.d/w);
	return part1*part2*part3;
}

double Func::fgs() {
	return pow((n-2.*d)/(n-2.*w*r), 2.);
}

double Func::fdgsdm() {
	return (-2.*(-2.*(-m + n/2.) + 2.*(m + n/2.))*Power(-2.*d + n,2))/Power(-2.*(-m + n/2.)*(m + n/2.) + n,3);
}

double Func::fdgsdd() {
	return (-4.*(-2.*d + n))/Power(-2.*(-m + n/2.)*(m + n/2.) + n,2);
}

double Func::fdgsdn() {
	return (-2.*(1. - n)*Power(-2.*d + n,2))/Power(-2.*(-m + n/2.)*(m + n/2.) + n,3) +
			   (2.*(-2.*d + n))/Power(-2.*(-m + n/2.)*(m + n/2.) + n,2);
}

double Func::fdgtdm() {
	return -(((Sqrt(((1. + d - n)*(1. + m - n/2.))/(1. - m - n/2.)) + Sqrt((d*(-m + n/2.))/(m + n/2.)))*
	        (-2*(-m + n/2.) + 2*(m + n/2.))*(Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) +
	          Sqrt((d*(m + n/2.))/(-m + n/2.)))*(-2*d + n))/Power(-2*(-m + n/2.)*(m + n/2.) + n,2)) +
	   ((Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) + Sqrt((d*(-m + n/2.))/(m + n/2.)))*
	      ((-(((1 + d - n)*(1 - m - n/2.))/Power(1 + m - n/2.,2)) - (1 + d - n)/(1 + m - n/2.))/
	         (2.*Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.))) +
	        (d/(-m + n/2.) + (d*(m + n/2.))/Power(-m + n/2.,2))/(2.*Sqrt((d*(m + n/2.))/(-m + n/2.))))*(-2*d + n))/
	    (-2*(-m + n/2.)*(m + n/2.) + n) + ((((1 + d - n)/(1 - m - n/2.) +
	           ((1 + d - n)*(1 + m - n/2.))/Power(1 - m - n/2.,2))/
	         (2.*Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.))) +
	        (-((d*(-m + n/2.))/Power(m + n/2.,2)) - d/(m + n/2.))/(2.*Sqrt((d*(-m + n/2.))/(m + n/2.))))*
	      (Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) + Sqrt((d*(m + n/2.))/(-m + n/2.)))*(-2*d + n))/
	    (-2*(-m + n/2.)*(m + n/2.) + n);
}

double Func::fdgtdd() {
	return (-2*(Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) + Sqrt((d*(-m + n/2.))/(m + n/2.)))*
		      (Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) + Sqrt((d*(m + n/2.))/(-m + n/2.))))/
		    (-2*(-m + n/2.)*(m + n/2.) + n) + ((Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) +
		        Sqrt((d*(-m + n/2.))/(m + n/2.)))*((1 - m - n/2.)/
		         (2.*Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.))*(1 + m - n/2.)) +
		        (m + n/2.)/(2.*(-m + n/2.)*Sqrt((d*(m + n/2.))/(-m + n/2.))))*(-2*d + n))/
		    (-2*(-m + n/2.)*(m + n/2.) + n) + (((1 + m - n/2.)/
		         (2.*(1 - m - n/2.)*Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.))) +
		        (-m + n/2.)/(2.*Sqrt((d*(-m + n/2.))/(m + n/2.))*(m + n/2.)))*
		      (Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) + Sqrt((d*(m + n/2.))/(-m + n/2.)))*(-2*d + n))/
		    (-2*(-m + n/2.)*(m + n/2.) + n);
}

double Func::fdgtdn() {
	return -(((Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) + Sqrt((d*(-m + n/2.))/(m + n/2.)))*
	        (Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) + Sqrt((d*(m + n/2.))/(-m + n/2.)))*(1 - n)*
	        (-2*d + n))/Power(-2*(-m + n/2.)*(m + n/2.) + n,2)) +
	   ((Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) + Sqrt((d*(-m + n/2.))/(m + n/2.)))*
	      (Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) + Sqrt((d*(m + n/2.))/(-m + n/2.))))/
	    (-2*(-m + n/2.)*(m + n/2.) + n) + ((Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) +
	        Sqrt((d*(-m + n/2.))/(m + n/2.)))*((((1 + d - n)*(1 - m - n/2.))/(2.*Power(1 + m - n/2.,2)) -
	           (1 + d - n)/(2.*(1 + m - n/2.)) - (1 - m - n/2.)/(1 + m - n/2.))/
	         (2.*Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.))) +
	        (d/(2.*(-m + n/2.)) - (d*(m + n/2.))/(2.*Power(-m + n/2.,2)))/(2.*Sqrt((d*(m + n/2.))/(-m + n/2.))))*
	      (-2*d + n))/(-2*(-m + n/2.)*(m + n/2.) + n) +
	   (((-(1 + d - n)/(2.*(1 - m - n/2.)) + ((1 + d - n)*(1 + m - n/2.))/(2.*Power(1 - m - n/2.,2)) -
	           (1 + m - n/2.)/(1 - m - n/2.))/(2.*Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.))) +
	        (-(d*(-m + n/2.))/(2.*Power(m + n/2.,2)) + d/(2.*(m + n/2.)))/(2.*Sqrt((d*(-m + n/2.))/(m + n/2.))))*
	      (Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) + Sqrt((d*(m + n/2.))/(-m + n/2.)))*(-2*d + n))/
	    (-2*(-m + n/2.)*(m + n/2.) + n);
}

/*
 *********************************************************
 ******************* CLASS Equations *********************
 *********************************************************
*/

class Equa {
private:
	parameters p; // parameters
	variables v; // variables
	Func *f;

	void onePoint(variables &y, int i, int j);

public:
	Equa(parameters &p, variables &v) : p(p), v(v) { f = new Func(p, v); }
	~Equa() { delete f; }

	void setPar(parameters &p) {
		this->p = p;
		f->set(p);
	}

	void run(variables &w);
};

void Equa::onePoint(variables &w, int i, int j) {
	double x = (double)i*2.*M_PI/p.L;
	double y = (double)j*2.*M_PI/p.L;

	f->refresh(x, y);
	if(DEEP_DEBUG) { f->printLambda(); }

	double Etan[4] = { (1-tanh(p.beta*f->E[0]/2.))/2., (1-tanh(p.beta*f->E[1]/2.))/2., (1-tanh(p.beta*f->E[2]/2.))/2., (1-tanh(p.beta*f->E[3]/2.))/2. };

	for(int i=0; i<4; i++) {
		w.chiAB += Etan[i]*f->dEdlchiAB[i];
		w.chiS += Etan[i]*f->dEdlchiS[i];
		w.chiT += Etan[i]*f->dEdlchiT[i];
		w.DeltaS += Etan[i]*f->dEdlDeltaS[i];
		w.DeltaT += Etan[i]*f->dEdlDeltaT[i];
		w.mu += Etan[i]*f->dEdln[i];
		w.mAF += Etan[i]*f->dEdlmAF[i];
		// w.d w tym miejscu jeszcze nie
	}
}

void Equa::run(variables &w)
{
	clear(w);
	variables temp; clear(temp); // temporary vector using during summation
	unsigned int max = p.L/2; // attention: L should be even!

	// SERIAL CASE
	// SMART WAY OF SUMMATION - I use symmetry k_x,y <-> -k_x,y
	// * edges (2 times)
	for(unsigned int i=1; i<max; i++) {
		onePoint(temp, i, 0);
		onePoint(temp, 0, i);
		onePoint(temp, i, max-i);
	}
   addTo(temp, w, 2);
   clear(temp);
   // * interior area (4 times) - SERIAL CASE
   for(unsigned int i=1; i<max; i++) {
		for(unsigned int j=1; j<i; j++) {
			onePoint(temp, max-i, j);
     	}
	}
	addTo(temp, w, 4);
	clear(temp);
	// * missing corners
	onePoint(w, max, 0);
	onePoint(w, 0, max);
	onePoint(w, 0, 0);

	// what was out of sumation
	double L2 = (double)(p.L*p.L);
	w.chiAB /= L2;
	w.chiS /= L2;
	w.chiT /= L2;
	w.DeltaS /= L2;
	w.DeltaT /= L2;
	w.mu /= L2;
	w.mAF /= L2;

	w.chiAB += 8.*v.chiAB;
	w.chiS += 8.*v.chiS;
	w.chiT += 8.*v.chiT;
	w.DeltaS += 8.*v.DeltaS;
	w.DeltaT += 8.*v.DeltaT;
	w.mu += 0.5*p.n;
	w.mAF += 0.5*v.mAF;

	// calculate dd
	w.d = 8.*(p.t*v.chiAB + p.tPrim*v.chiS)*f->dgtdd + p.J*(-3.*v.chiAB*v.chiAB + v.DeltaS*v.DeltaS \
				+ 3.*v.DeltaT*v.DeltaT + v.mAF*v.mAF/2.)*f->dgsdd + 2.*v.d*p.U;

	//if(w.chiAB != w.chiAB) { w.chiAB = 1000; }
	//if(w.chiS != w.chiS) { w.chiS = 1000; }
	//if(w.chiT != w.chiT) { w.chiT = 1000; }
	//if(w.DeltaS != w.DeltaS) { w.DeltaS = 1000; }
	//if(w.DeltaT != w.DeltaT) { w.DeltaT = 1000; }
	//if(w.mu != w.mu) { w.mu = 1000; }
	//if(w.mAF != w.mAF) { w.mAF = 1000; }
}

/////////////////////////////////////////////////////////////////////////////////////////////////



/*
 *********************************************************
 ************************* GSL ***************************
 *********************************************************
*/

	// parameters  p.n  p.J  p.t  p.tPrim  p.U  p.beta  p.L
	// variables  v.chiAB  v.chiS  v.chiT  v.DeltaS  v.DeltaT  v.mu  v.mAF  v.d

int set (const gsl_vector * x, void *params, gsl_vector * f) {
	parameters p = *((struct parameters *) params);

	variables v;
	v.chiAB = gsl_vector_get (x, 0);
	v.chiS = gsl_vector_get (x, 1);
	v.chiT = gsl_vector_get (x, 2);
	v.DeltaS = gsl_vector_get (x, 3);
	v.DeltaT = gsl_vector_get (x, 4);
	v.mu = gsl_vector_get (x, 5);
	v.mAF = gsl_vector_get (x, 6);
	v.d = gsl_vector_get (x, 7);

	if(DEEP_DEBUG) {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15ee\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3), gsl_vector_get (x, 4), gsl_vector_get (x, 5),
				gsl_vector_get (x, 6), gsl_vector_get (x, 7));
	}

	Equa E(p, v);
	variables w;
	E.run(w);

	gsl_vector_set (f, 0, w.chiAB);
	gsl_vector_set (f, 1, w.chiS);
	gsl_vector_set (f, 2, w.chiT);
	gsl_vector_set (f, 3, w.DeltaS);
	gsl_vector_set (f, 4, w.DeltaT);
	gsl_vector_set (f, 5, w.mu);
	gsl_vector_set (f, 6, w.mAF);
	gsl_vector_set (f, 7, w.d);

	if(DEEP_DEBUG) {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3), gsl_vector_get (x, 4), gsl_vector_get (x, 5),
				gsl_vector_get (x, 6), gsl_vector_get (x, 7));
		printf("###  %24.15e# %24.15e# %24.15e# %24.15e# %24.15e# %24.15e# %24.15e# %24.15e\n",
			gsl_vector_get (f, 0), gsl_vector_get (f, 1), gsl_vector_get (f, 2),
			gsl_vector_get (f, 3), gsl_vector_get (f, 4), gsl_vector_get (f, 5),
			gsl_vector_get (f, 6), gsl_vector_get (f, 7));
		printf("\n");
	}

	return GSL_SUCCESS;
}

// parameters  p.n  p.J  p.t  p.tPrim  p.U  p.beta  p.L
void print_label(parameters &p) {
	printf("\nn=%f, J=%f, t=%f, tPrim=%f, U=%f, beta=%f, L=%u\n", p.n, p.J, p.t, p.tPrim, p.U, p.beta, p.L);
	printf("iter \tchiAB \t\tchiS \t\tchiT \t\tDeltaS \t\tDeltaT \t\tmu \t\tmAF \t\td\n");
}

void print_state (size_t iter, gsl_multiroot_fsolver * s)
{
	printf("%3u %10.3e (%2d) %10.3e (%2d) %10.3e (%2d) %10.3e (%2d) %10.3e (%2d) %10.3e (%2d) %10.3e (%2d) %10.3e (%2d)\n",
			(unsigned int)iter, gsl_vector_get (s->x, 0), (int)(log(abs(gsl_vector_get (s->f, 0)))/log(10)),
			gsl_vector_get (s->x, 1), (int)(log(abs(gsl_vector_get (s->f, 1)))/log(10)),
			gsl_vector_get (s->x, 2), (int)(log(abs(gsl_vector_get (s->f, 2)))/log(10)),
			gsl_vector_get (s->x, 3), (int)(log(abs(gsl_vector_get (s->f, 3)))/log(10)),
			gsl_vector_get (s->x, 4), (int)(log(abs(gsl_vector_get (s->f, 4)))/log(10)),
			gsl_vector_get (s->x, 5), (int)(log(abs(gsl_vector_get (s->f, 5)))/log(10)),
			gsl_vector_get (s->x, 6), (int)(log(abs(gsl_vector_get (s->f, 6)))/log(10)),
			gsl_vector_get (s->x, 7), (int)(log(abs(gsl_vector_get (s->f, 7)))/log(10)));
}

// variables  v.chiAB  v.chiS  v.chiT  v.DeltaS  v.DeltaT  v.mu  v.mAF  v.d
void print(FILE * stream, double x, double y, variables &v) {
	fprintf(stream, "%f %f %.15e %.15e %.15e %.15e %.15e %.15e %.15e %.15e       ",
			(float)x, (float)y, v.chiAB, v.chiS, v.chiT, v.DeltaS, v.DeltaT, v.mu, v.mAF, v.d);
	//fprintf(stream, "%.15e %.15e %.15e %.15e %.15e %.15e %.15e\n",
	//		v.lchiAB, v.lchiS, v.lchiT, v.lDeltaS, v.lDeltaT, v.ln, v.lmAF);
	fflush (stream);
}

int aPoint(parameters &p, variables &v) {
	const gsl_multiroot_fsolver_type *T;
	gsl_multiroot_fsolver *s;
	int status;
	size_t iter = 0;
	const size_t n = NO_OF_EQUATIONS;
	gsl_multiroot_function f = {&set, n, &p};

	double x_init[n];
	x_init[0] = v.chiAB;
	x_init[1] = v.chiS;
	x_init[2] = v.chiT;
	x_init[3] = v.DeltaS;
	x_init[4] = v.DeltaT;
	x_init[5] = v.mu;
	x_init[6] = v.mAF;
	x_init[7] = v.d;

	gsl_vector *x = gsl_vector_alloc (n);
	for(size_t i=0; i<n; i++) { gsl_vector_set (x, i, x_init[i]); }

	T = gsl_multiroot_fsolver_hybrids;
	s = gsl_multiroot_fsolver_alloc (T, n);
	gsl_multiroot_fsolver_set (s, &f, x);

	if(DEBUG) {
		print_label(p);
		print_state (iter, s);
	}

	do {
		iter++;
		status = gsl_multiroot_fsolver_iterate (s);
		if(DEBUG) { print_state (iter, s); }
	    if(status) break;	// check if solver is stuck
	    status = gsl_multiroot_test_residual (s->f, PRECISION);
	}
	while (status == GSL_CONTINUE && iter < MAX_ITERATION);

	v.chiAB = gsl_vector_get (s->x, 0);
	v.chiS = gsl_vector_get (s->x, 1);
	v.chiT = gsl_vector_get (s->x, 2);
	v.DeltaS = gsl_vector_get (s->x, 3);
	v.DeltaT = gsl_vector_get (s->x, 4);
	v.mu = gsl_vector_get (s->x, 5);
	v.mAF = gsl_vector_get (s->x, 6);
	v.d = gsl_vector_get (s->x, 7);

	if(DEBUG) printf ("(aPoint) status = %s", gsl_strerror(status));
	gsl_multiroot_fsolver_free (s);
	gsl_vector_free (x);
	return status;
}

/*
 *********************************************************
 ******************** INTERPOLATION **********************
 *********************************************************
*/

//double extrapolation(double &x, double &y, double point) {
//	gsl_interp_accel *acc = gsl_interp_accel_alloc ();
//	gsl_spline *spline = gsl_spline_alloc (gsl_interp_cspline, POINTS_TO_INTERP);
//
//	gsl_spline_init (spline, x, y, POINTS_TO_INTERP);
//	double w = gsl_spline_eval (spline, point, acc);
//
//	gsl_spline_free (spline);
//	gsl_interp_accel_free (acc);
//	return w;
//}

void shake(variables &v, double a) {
	v.chiAB += a*((double)rand()/(RAND_MAX+1.0));
	v.chiS += a*((double)rand()/(RAND_MAX+1.0));
	v.chiT += a*((double)rand()/(RAND_MAX+1.0));
	v.DeltaS += a*((double)rand()/(RAND_MAX+1.0));
	v.DeltaT += a*((double)rand()/(RAND_MAX+1.0));
	v.mu += a*((double)rand()/(RAND_MAX+1.0));
	v.mAF += a*((double)rand()/(RAND_MAX+1.0));
	v.d += a*((double)rand()/(RAND_MAX+1.0));
}

void interpolation1(variables &v, variables &old) {
	v.chiAB += v.chiAB - old.chiAB;
	v.chiS += v.chiS - old.chiS;
	v.chiT += v.chiT - old.chiT;
	v.DeltaS += v.DeltaS - old.DeltaS;
	v.DeltaT += v.DeltaT - old.DeltaT;
	v.mu += v.mu - old.mu;
	v.mAF += v.mAF - old.mAF;
	v.d += v.d - old.d;
}

/*
 *********************************************************
 ************************ MAIN ***************************
 *********************************************************
*/

#define POINTS_TO_INTERP 5
#define START 20
#define STOP 0
#define OUTPUT_FILE "test.dat"
#define frand() ((double)rand()/(RAND_MAX+1.0))

#define U_MIN 2.
#define U_MAX 20.
#define U_RESOLUTION 36.
#define N_MIN 0.65
#define N_MAX 0.97
#define N_RESOLUTION 64.

// it is what we are running really!
int yttSGA (void) {
	srand (time(NULL));

	// parameters  p.n  p.J  p.t  p.tPrim  p.U  p.beta  p.L
	struct parameters p;
	p.n = 0.93;
	p.J = 1./3.;
	p.t = -1.;
	p.tPrim = 0.25;
	p.U = 10;
	p.beta = 500;
	p.L = 512;

	// variables  v.chiAB  v.chiS  v.chiT  v.DeltaS  v.DeltaT  v.mu  v.mAF  v.d
	struct variables v;
	v.chiAB = -2.026352022056961e-01;
	v.chiS = 3.836555826800296e-03;
	v.chiT = -8.661321041808190e-11;
	v.DeltaS = -5.172879501188971e-13;
	v.DeltaT = 6.788427458416529e-12;
	v.mu = 2.136173396478271e-01;
	v.mAF = 2.739235294447910e-10;
	v.d = 2.294560104496349e-01;

	//SPRAWDZIĆ (dla t'=0)
	//20.000000 0.990000 -2.026352022056961e-01 3.836555826800296e-03 -8.661321041808190e-11 -5.172879501188971e-13 6.788427458416529e-12 2.136173396478271e-01 2.739235294447910e-10 2.294560104496349e-01       	1.099468846228582e+00
	//20.000000 0.990000 -2.021842806537678e-01 3.699863835479314e-03 2.605094968534017e-02 -2.805624678085749e-12 3.646263684048568e-12 1.555859324865105e-01 -1.204080044089004e-01 2.305526366245958e-01       	1.115681773074186e+00


	FILE * pFile;
	pFile = fopen(OUTPUT_FILE, "w");
	fprintf(pFile, "# n=%f, J=%f, t=%f, tPrim=%f, U=%f, beta=%f, L=%u\n", p.n, p.J, p.t, p.tPrim, p.U, p.beta, p.L);
	//fprintf(pFile, "zmiana n od 0 do 0.35\n");
	fprintf(pFile, "# x \t\t\tchiAB \t\t\tchiS \t\t\tchiT \t\t\tDeltaS \t\t\tDeltaT \t\t\tmu \t\t\tmAF \t\t\t  d");


	// sprawdzam jeden punkt
//	aPoint(p, v);
//	return 0;


	variables copy1;
	copyTo(copy1, v);
	Func *fun = new Func(p,v);

	// sprawdzam jeden punkt (ciągle)
	while(1) {

		printf("\nJedziemy\n");
		int status = aPoint(p, v);
		if(status == GSL_SUCCESS) {
			fprintf(pFile, "\n");
			print(pFile, p.U, p.n, v);
			// energia dla znalezionych punktow
            fun->set(p, v);
            double energy = fun->Energy();
            fprintf(pFile, "\t%20.15e", energy);
            printf(" E = %20.15e\n", energy);
			fflush(pFile);
		}

		copyTo(v, copy1);
		//shake(v, 0.2);
		v.chiAB = 2.*frand()-1.;
		v.chiS = 2.*frand()-1.;
		v.chiT = 2.*frand()-1.;
		v.DeltaS = 2.*frand()-1.;
		v.DeltaT = 2.*frand()-1.;
		v.mu = 2.*frand()-1.;
		v.mAF = 2.*frand()-1.;
		v.d = frand();
	}

	fclose (pFile);
	delete fun;
	return 0;
}
/*
	Func *fun = new Func(p,v);

	p.U = U_MAX;
	p.n = N_MAX;
	variables old, old2;
	aPoint(p, v);
	copyTo(old, v);

//	FILE *gnuplot_pipe=popen("gnuplot", "w");
//	fprintf(gnuplot_pipe, "load 'init.gnu'\n");
//	fflush(gnuplot_pipe);

	double u_step = (U_MAX - U_MIN)/U_RESOLUTION;
	double n_step = (N_MAX - N_MIN)/N_RESOLUTION;

	for(double j=U_MAX; j>=U_MIN-1e-5; j-=u_step/2.) {
		p.n = N_MAX;
		p.U = j;
		copyTo(v, old);
		aPoint(p, v);

		interpolation1(v, old);
		shake(v, 1e-5);

		copyTo(old, v);
		copyTo(old2, v);
		for(double i=N_MAX; i>=N_MIN-1e-5; i-=n_step) {
			p.n = i;

			int status = aPoint(p, v); // aPoint take initiate point from v vector and put results back to v
			printf("%f %f : satus = %s\n", j, i, gsl_strerror(status));
			if(status != GSL_SUCCESS) {
				copyTo(v, old2);
            	shake(v, 1e-3);
				continue;
			}
			else {
				fprintf(pFile, "\n");
				print(pFile, j, i, v);

				// energia dla znalezionych punktow
            	fun->set(p, v);
            	double energy = fun->Energy();
            	fprintf(pFile, "\t%20.15e", energy);
            	printf(" E = %20.15e", energy);

				fflush(pFile);

				variables temp;
				copyTo(temp, v);

				interpolation1(v, old2);
				shake(v, 1e-5);

				copyTo(old2, temp);
			}
			// aPoint powinno zwracać status. Jeśli nie tak to np. cofnąć się? (zagęścić pomiary)
			// dodać opcję ekstrapolacji
		}
			fprintf(pFile, "\n");
			fflush(pFile);
//			fprintf(gnuplot_pipe, "load 'test2D-tJU-vol.gnu'\n");
//			fflush(gnuplot_pipe);

	}
	return 0;
}

*/
